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Abstract. We study the application of various forms of the coupled cluster method 
to systems with paired fermions. The novel element of the analysis is the study of 
the breaking and eventual restoration of particle number in the CCM variants. We 
specifically include Arponen's extended coupled cluster method, which describes the 
normal Hartree-Fock-Bogoliubov mean field at lowest level of truncation. We show 
that all methods converge to the exact results as we increase the order of truncation, 
but that the breaking of particle number at an intermediate level means that this 
convergence occurs in a surprising way. We argue that the most straightforward form 
of the method seems to be the most stable approach to implement for realistic (large 
number of particles) pairing problems 
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1. Introduction 

The Coupled Cluster Method (CCM) has a long and venerable history in both Physics 
and Chemistry. It is a method that can be designed to incorporate the normal mean-field 
approximation, but in principle it also allows one to make systematic improvements, 
and in principle converges to the exact solution. After initial work by Coester and 
Kmmel PQ, and Kmmel and his group [2] with important applications in nuclear 
physics, this was picked up in quantum chemistry by Cizek and Paldus [3j HJ 13 E] 
and took off from there (see Ref. [7J for a recent review). About 10 years ago, the 
method made a renewed appearance in nuclear physics in the work of Mihaila and 
Heisenberg [8j HJ [TOj [11], and in a more advanced application in the work of Dean 
et al ^T^MUMMTilMUMM)lEIl&M- Atomic nuclei are indeed interesting 
and complicated, even more so once we start thinking about heavier open shell nuclei, 
since we then have to take into account nuclear superfluidity, described by pairing of 
nucleons of opposite spins, as well as the "normal" problems of nuclear physics (short 
range repulsion, intermediate attraction, non-central forces). The normal way to deal 
with such nuclei using involves mean-field techniques, but one would like to use a fully 
microscopic method such as the CCM to tackle the full many-body problem, including 
the pairing part. Since we have some understanding of the nuclear physics complications, 
in this paper we shall concentrate on the pairing problem, to see how the CCM can be 
applied effectively to this aspect of the nuclear landscape. 

Of course pairing is also extremely important in atomic (fermionic) condensates, 
where we have the added advantages that experimentalists can tune the s-wave 
scattering length, which determines the interaction in the pairing channel, to almost 
any value required [23]. This has allowed one to study the BEC to BCS transition, 
and has given many-body theory access to a whole new set of data that require an 
explanation. The extreme elegance and control of the experiments means that the data 
requires high accuracy calculations, and is also able to probe collective modes, some 
of which are beyond simple mean-field plus harmonic fluctuations (RPA) calculations. 
One technique that gives one easy access to such modes is again the CCM. This gives a 
second impetus to this study. 

Little work has been done within CCM for problems with pairing forces. The so- 
called normal CCM has been used by Emrich and Zabolitzky fM\ 123] . but this work 
seems to have been largely forgotten, and it is probably not as general as required. If we 
are interested in systems with strong correlations, we may want to use methods that go 
beyond simple mean-field theory. The coupled cluster method has shown to be a very 
powerful technique to make such improvements, but we have little intuition what the 
best way to apply the method is. The work of Arponen [2S] and Arponen and Bishop 
[271 EH ESI EDI EU [32] provides one way; the use of Br ueckner /maximum overlap orbitals 
[33] . where we start from the BCS wave function, and express all corrections in terms 
of BCS quasiparticles, is an alternative. 

In this paper we shall investigate both these methods for a simplified model, with 
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an emphasis on particle number-so we concentrate on a single one of the aspects that 
are important for finite nuclei. Before analyzing this model, let us first describe the 
methods used succinctly. 

2. The normal CCM 

Let us look at the traditional (normal) Coupled Cluster Method (NCCM) first, both to 
make contact with other notations used in the field, and for later reference. 

We shall use the standard "Bochum/Manchester notation", where S denotes the 
cluster operator (called "T" in most chemically inspired literature). We also use the 
"SUB(n)" notation, where SUB(l) describes the truncation of the cluster operator to a 
single quasiparticle excitation, SUB (2) single plus double (SD), etc. Following Arponen 
[27?] . we start from a functional, 

0(s,s) = (0|O|0) = (0 o |(l + 5)e- 5 Oe 5 |0o), (1) 

which is used to evaluate the expectation value of any operator including the 
Hamiltonian. This has interesting linking properties: all S operators contract with 
O, and the S has to link with at least one S or the O operator. 

The states used in the expectation value are not Hermitian conjugates, but are 
intermediate normalised, 

(0o|(l + 5)e-V|0 o ) = l- (2) 
The resulting expectation value is a finite polynomial once we impose the condition 

S = J2 SlC i> S = J2~^Ci, (3) 
/ i 

with the generalised annihilation operators 

Ci |0o) = 0. (4) 

The reason for the polynomial nature of the result is that the nested commutator 
expansion 

e~ s Oe s = + [0,S] + ^[[0,S},S] + ... 

terminates at finite order if the operator O contains a finite number of operators. 

Normal CCM is usually written in a slightly different form; if we extremise the 
energy (as expectation value of the Hamiltonian) with respect to Si and sj, we only 
have to solve for sf. 

(MCie-^He 8 ^) = 0, (5) 
<0o|e- s °W c< U> =E, (6) 

where we use "eq" to denote the equilibrium solution of Of course to calculate the 
expectation value of any other operator we require as well, which can be found from 
a set of inhomogeneous linear equations 

(0o|(l + ^) e - 5cq [if,C / t ]e 5cq |0o) = O. 
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Finally, by expanding the action 

A[s, s} = J dt{4\{\ + S)e~ s (H - id t )e s \<f>), (7) 

to second order about the equilibrium we get the "Harmonic approximation" , also called 
EOM-CC [M] , with classical action, derived by shifting sj — > s e j q + sj 

A = J {i^h + s^h + ^isj [d Sl d Sj H(s,r q )\s=s^] 

+ S!Sj [d Sl d SJ H(s, s)\ s=s ^] \dt. (8) 



(The linear term, a total derivative, is a boundary term and can thus normally 
be removed from the action.) Varying with respect to £/(£), assuming a harmonic 
dependence of sj, Sj = eie tujt , we get 

[d Sl d Sj H(s,s eq )\ s=s c q ]ej = ua . (9) 

If we wish to get the tilde eigenstates as well, we need to vary with respect to sj(t), 
using sj = fie lU)t 

[d Sl d Sj H(s, S eq )| s= , cq ] ej + [d Sj d Sl H(s, s )| s=s e q ] fj = -ufr. (10) 

We can rewrite this as the eigenvectors of what looks like an RPA excited state problem 
[35J (which is due to the fact that this is the diagonalisation of a quadratic Hamiltonian) 

d Sl d S} H(s,s C(l )\ s=s ^ d Sj d Sl H(s,s)\ s=s c Cl \ = fo -I \ 
d Sj d~ Sl H(s,s)\ s=sC , J \ I J ' 1 1 

where we have used the "double height column vector" 

*=(l}- (12) 




For this particular matrix the eigenvalues are completely determined by the 
diagonalisation of the lower left block d Sj dsjH(s, s )| s=s eq, which is the matrix of 
derivatives of the traditional CCM equations, which are given by the Eqs. (jSJ), and 
are used to determine the coefficients s/. 



3. The extended CCM 



The extended CCM starts from a similar functional as Eq. ([I]) , but now based on a 
double-exponential form of the bra states, 

0(5, s) = (4>\0\4>) = (0 o |eV 5 Oe 5 |0o). (13) 

The operators remain unchanged, i.e., they are defined in Eqs. ([3]) and (j3J). Unlike 
the ordinary NCCM, it has not been widely applied since it is much more difficult to 
evaluate. Here we discuss a quasiparticle approach that seems to hold some promise. 
We shall in essence make intimate use of a link to the theory of generalised coherent 
states, and the fact that Eq. (|T3"|) is identical to the Dyson-Maleev boson mapping of 
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Lie algebras [36] when we can assign a Lie algebra to the operators Ci,Cj and their 
commutators. 

When we minimize the energy, we can now no longer separate the steps in the 
determination of sj and sj, but have to combine the equations 

(0o|C7e* c V 5O W c >o> =0, 
(0o|e^ q e- scq [ J ff,C)]e 5C > o ) = O ) 

(Me^e^Hen^o) = E. (14) 

Once we have solved these, we can look for the excited states by solving the harmonic 
problem, sj — >■ + s/, etc., 



A= j{ 1 X>o|e%j|0o>s/ + l -s lSj [d Sl d Sj H(s, ~s c 

+ s lSj [d Sl d SJ H(s,s e(l )\ s=s e^ dt.} (15) 

Varying with respect to Sj(t), assuming a harmonic dependence of sj, sj = eie luJt , we 
get 



d Sl d Sj H(s,s) d Sj d Sl H(s,s 
d~ Sl d Sj H(s,s) d Sl d Sj H(s,s 



d Sl (<p \e § C\\<t> ) 




(16) 



s=i a ...Lt,,, — , > (17) 



cq 



where the symplectic matrix S can be brought to canonical form by replacing the set 
of amplitudes {5j} by the new set {cj} which solves the nonlinear equations 

<7j = <0o|e 5 C}|<M, (18) 

and expanding the energy in these new amplitudes. This is not necessary to solve the 
RPA, however, and for certain mixed calculations (see below) such a transformation is 
not easily done, and we may want to choose to work with the form given in (I16|17p . 



4. ECCM SUB(l) approximation and HFB 

One of the natural ways to generalise the Hartree-Fock method to problems with pairing, 
and the BCS approach to problems with particle-hole interactions is the Hartree-Fock- 
Bogoliubov method. We still use a single Slater-determinant wave function, but the state 
has no definite particle number, like the BCS state. As shown below, we can capture 
this information in a generalised density matrix [35J. There is a CCM approximation 
that is equivalent to this state: 

The ECCM SUB(l) approximation for a finite fermionic system, where the reference 
state is a single Slater determinant of occupied orbitals h 

i0o)=n a ii°) ( i9 ) 

h 
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is denned by the one-body S operatoJ|| 

S = s vh a\a h + sPP ' a H' + E skh ' a h' a h, (20) 

ph pp' hh' 

where h labels the single-particle states occupied in the reference state, and p the 
unoccupied ones. (In other words, this is the "singles" truncation of quantum chemistry, 
adapted to the pairing problem where we no longer preserve particle number). We now 
evaluate ({TBI by first using the nested commutator expansion for e~ s He s . In calculating 
the commutator of S with the Hamiltonian, either one or both of the operators in 
S contract with H. If both contract, no further contractions are possible. Having 
performed all possible contractions with S we need to now calculate the final contractions 
with S. In ECCM this will either tie together two S's, one S and link the other operator 
to H, or link both its indices to H. We thus get a factorisation of S and S into general 
objects with two external indices, essentially one-body densities, which is exactly the 
Hartree-Fock-Bogoliubov (HFB) factorisation, if we identify the linked strings of S and 
5"s with the normal and abnormal densities. This can of course be made explicit: the 
algebra is rather lengthy, but see below for a simplified (BCS) example. 



5. Simplifying SUB(l) 



Let us look at the application of the extended CCM to a pure pairing problem (without 
particle-hole excitations) in the SUB1 approximation. We start from the standard 
definitions 

Sp Sh 



S— s P a vt a -vi + / ] s hO'hjO'-ht, 
p h 



(21) 
(22) 



We also write 



e S |0o), 



— 

s h 



i> |e 5 e s 



(23) 



We now analyse the normal and abnormal densities in the standard way. 



The normal density pi a j a ' = 



Pia,j<r' 3aa' 



equals 



SpSp 



(24) 



h p 

and the abnormal densities require a bit more work, but after some algebra can be 
expressed as 



X Thus the SUB(n) calculation has an n-body S operator. 
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| Q'ja'Q'ia | 



[Sat S <r'i - Scri,Sa't\ { h^ij + S p (l - S p S p )5^ 



(25) 



The bar denotes the time-reversed state 



5.1. Generalised density 

From the standard block-matrix definition of the generalised density, 



p K 
-K* I - p* 



(26) 



we can easily check that 1Z 2 = 1Z. Expanding out this projector condition, we get the 
four matrix conditions 

2 * 
p — KK = p, 

pK — K,p* = 0, 

P K — K p = (J, 

p — K K = p . (27) 

Let us show explicitly how one of these works for the forms found above. First notice 
that 



Thus 



£ ~ 3hS ^ ~ ShSh + E 5P ij(~ s v s p) C 1 - 5 p s 



^■(SpSp) 



— = o. 



(29) 



h p 
which add up to p. 

Of course, the other three elements can be evaluated in a similar way. This shows 
that we have a general non-Hermitian classical mapping of the general density matrix. 
[Discuss mapping issues] 

5.2. Canonical form 

In the literature (e.g., Ref. [35]) we can find the canonical form for the density matrix 



as 



p = ^<W + u\8 hth! 

K = K* = U k V k 8 k -y. 



(30) 
(31) 



Our expressions should be of the same form, since we have assumed the diagonal form 
of S. Actually, we pay a price here for unnecessarily starting with a Hartree-Fock state; 
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all calculations above are valid with respect to the vacuum as well, where we recover a 
simpler form of BCS (which has maximal asymmetry of k\). From the relations 



P = v l y k ,k' + %f) ' ( 32 ) 

K = K * = UkVk S k —^ (33) 

where we have rewritten the basis as 

k = k |, k = -k I, (34) 

we thus conclude that 

s k = v k u k = -sm2tp k , (35) 

Sk = v k /u k = cot^. (36) 



The energy, which is equal to the BCS energy, and is thus fully variational, is 

I(s, s) = 2e k s k s k + v klk2klk2 s kl s kl s k2 s k2 + v k2k2klkl s kl (1 — s kl s kl ) s k2 . (37) 

This needs to be extended with a Lagrange multiplier for particle number-since we 
break particle number we need to make sure that at least the average particle number 
is correct. At this level if truncation, we get a term —A (s k s k — Nq) in the functional. 
[Needs to be mentioned earlier as well] 



5.3. Quasiparticle basis 



It is interesting to note that the parametrisation above is naturally linked to a bi- 
canonical operatorial basis (c 1 ^ d\ but {dj,, q} = 5fc,z), which is obtained by calculating 
the action of a k on the ket |0) and a\ on the bra (<f>\ (the idea to use bi-canonical 
operators traces back at least as far as Ref. [37]), in each case subtracting the result 
from the initial operatorjjj 

ak + [S, a k ] 



Cfe 



<1i 


= a k 


+ [S, a~ k ] 


= a 


k + s kd k , 


4 


= a l 


+ 


5,4 


+ 


's, 


S,a[ 




4 


= 4 


+ 




+ 




5,4 





This can be inverted to give 
ak = (1 



SkSk)Ck + s k d^, 



[1 - s k s k )c^ - s k d[, 



,t 



- SfcCfc + 



4 



;i - s k s k )a k - s k a^, 



s k s k )aL + Sk a k . 



(38) 



(39) 



and 



ak \4>) = s k(J 
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Using these relations we can transform the original Hamiltonian (there is no sum over 
spin indices, since we have chosen a Hamiltonian with S'-wave pairing only): 

k,a fcifc 2 fc 3 fc 4 

If we use a time-reversal invariant potential, Vk^k^ = Vk 3 k 4 k x k 2 i an d we also assume 
spin-balance, v klk2 k 3 k A = ^fc 2 fcifc 4 fc 3 , we have 

F qp = R + H 20 d ^_ + H 20 c _ Ck + H UJ kCk + #11^ 

+ T . d£ (it 4 c fe4 + Hl\ xxdldLS-cu 

k 1 k 2 k 3 k 4 fei fc 2 k s K4 k 1 k 2 k 3 k 4 fci fc 2 fc 3 fc 4 

+ <W 4 <4 2 % c ^ ( 41 ) 

which is the non-Hermitian analogue of the standard quasiparticle Hamiltonian [35J. As 
the normal quasiparticle Hamiltonian, it can be used to simplify calculations. 

5. 4- Higher order terms and Brueckner orbitals 

Of course we can calculate higher order CCM contributions as well; see below for some 
explicit examples. In general, one of the main problems in applying the extended CCM 
in this way is the proliferation of terms ("diagrams") containing contributions- 
these describe the rotation of any single particle state, external or internal, whereas in 
NCCM we only rotate external states. One way to avoid this proliferation, is to use 
the quasiparticle Hamiltonian (I4ip defined above, and express the correlation through 
an operator expressed in the quasiparticles-or for that fact an even higher-order 
operator- 

^ (2) = s feife2fc 3 A44 1 4 2 4 3 4 4 ' ^ (2) = ^ik^^i^^s ^ ( 42 ) 

i.e. 



/(s< 2) ,S< 2) ,s,S) 



! + ^» \e- Sm H«°e^m. (43) 



where the alternative shows that we can have either ECCM or NCCM for the extra 
terms. This corresponds to the use of what is called "Brueckner orbitals" in Quantum 
Chemistry, as long as we optimize the quasiparticle operators and the S^ n \ n > 2, at 
the same time. Since the inclusion of higher order terms is our main target, we shall 
look specifically at these two methods. 
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6. Single shell pairing 

6.1. Model investigation 

In order to analyse the various approaches, we study a model. The simplest model 
exhibiting many of the features we are interested in, is the single-shell pairing model. 
This has a Hamiltonian of the form [32] 

H = -G{A j A- N/2). (44) 

In the context used here, we interpret it as describing Q states, each with spin up or 
dowrjjj]. In both cases we can divide the states into two equal size sets labeled by k and 
k, which are then paired, 

n 

k=l 

a n 

n = j2 a W + Y, a l a k- ( 45 ) 

k=l k=l 

These operators form the well-known SI(2) quasispin algebra [35] 
[A, A f ] = 2(^/2 - N/2) = 2J , 

[J , At] = -At, (46) 
[Jo, A] =A. 

From the algebra, we can easily derive the result that the non-degenerate ground state 
for even particle number N is the state A^ N ^ 2 \0), and the energy of this groundstate is 

^ = - G ( fi -f)y- < 47 > 

6.2. COM 

In order to find the coupled-cluster approximation to the exact solution, we need to find 
the extremum of the constrained energy 

E(s, s) = (0\e- § e s [H - \(N - N )] e s \0). (48) 

We find to lowest order that s cannot depend on the magnetic quantum number, since 
there is no preferred direction, and thus 

E(s, s) = -A (2QsS - N ) - GQ(Q - l)ss(l - sns). 

Isolating the chemical potential term, we see that 

2n S § = (N) = N 0} (49) 

and thus 

j| In the original context it has states with angular momentum j = f2 — 1/2, and states of opposite 
angular momentum projection quantum numbers form pairs 
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We see a generic feature of the problem emerging: since we cast a variational calculation 
as a bi-variational one, introducing superfluous variables, there is no unique solution to 
s and s separately. The standard choice would be 

<A> = (At>, (51) 

where we find 



N 



20 -n' 



N (2n-N ) 

3 = V (2ny ' (52) 

but it is not necessary to make this choice. The SUB(l) approximation to the ground- 
state energy is thus 

E = . G ^^l(n-^], (53) 



2 o v 2 , 

which up to the correction of relative order I/O is the exact answer. 

6.2.1. U P article" EC CM We shall first study the ECCM based on the normal 
( "particle" ) operators, which, although normally a very lengthy calculation, is a simple 
calculation for this model, since we can use the quasispin algebra fj46|) to work out all 
results. In light of the fact that the exact solutions are of the form A^ 2 |0), we shall 
only consider the limited set of CCM states with structure 

|$)=exp^s n A+ n J |0>. (54) 

We have performed calculations with M up to 7. 

6.2.2. Quasiparticle ECCM The quasiparticle Hamiltonian is quite simple in this case. 
We use the quasiparticle number operator n = d\ck + dic%, and the quasiparticle pair 
operators 5 = c^Ck and 8^ = d^d^ to write 

H/G = Q /2 - (1 - ss)ss{VL - n) 2 

+ [(1 - s~sf + (ss) 2 ] [5^5 - (O - n)/2] 
+ sW + [S(l-s5)] 2 M 

- s(l - 2ss)5\n - n - 1) - 3(1 - s3)(l - 2s5)(Q -n- 1)5.(55) 

At the same time 

N = - (1 - 2ss)(n - n) +0 - 2s<5 f - 25(1 - ss)8, (56) 

and 

N 2 = tt 2 + 40s5 f + 405(1 - ss)5 

+ 8(1 - 88)88^6 

+ (1 - 2s5) 2 (0 - nf + 2(2s5(l - ss) - (1 - 2s5)0)(0 - n) 
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+ 4s W + 4 [S(l - s§)] 2 cW 

- 4s(l - 2ss)5\Q - n - 1) - 45(1 - sS)(l - 2ss)(n -n- 1)5.(57) 
We now calculate the expectation value of H — XN in the more general state 

(-\e §(M) e- s(M) (H-XN)e s(M) \-). (58) 
where we use the SUB(M) approximation, 

M M 

S (M) = J- i- s (") (5t) n , £W = ^ IgW^, (59) 

n=2 n ' n=2 n ' 

and 

|-) = e 5(1) |0), (-| = (0|e^ (1) e- 5(1) . (60) 
All relevant commutators can easily be done using the quasiparticle quasispin algebra, 

[5, (o" f ) m ] 
[n, (ft)™] 

[8n, (5T] 
[5\ (ft)™] 

Finally we need 

<-|eW 

6.3. Results 

6.3.1. ECCM based on particles We first wish to analyse the convergence of the ECCM 
based on the original operators. This has two reasons. First of all this will be a good 
test of the technology used in the implementation of the equations, which are all based 
on the algebraic implementation of the quasispin algebra-we shall see below why this is 
important. Also, we would like to understand what type of solutions we get, and how 
we approach the exact solution. 

With the help of Mathematica, we can obtain the full solution to the CCM 
equations, i.e., all real solutions to the coupled polynomial equations, see Fig. [TJ 
As we increase the truncation from the SUB(l) level-which is just the normal BCS 
approximation- to the SUB (2) level we see some new solutions appearing. The green 
line, which is not a very good solution in most cases, is actually exact at half filling. 
The blue line is a better approximation than BCS up to about half filling. The SUB (3) 
results give a rather busy picture, with lots of solutions, among which are exact solutions 
(for iV = 0, 2, 4, 6, 8). Finally for SUB(4), where ECCM coincides with NCCM, we find 
exactly four branches of the solution, but each physical solution lies on a different 



= (V) m -\tt-n-m-l), 
= 2m(5 t ) m , 

= m{5^) m -\Q - n -m+ l)(n + 2m) + 2m{5^) m 5, 
= m{m - l){5 ] ) m - 2 {Q - n - m + 2)(fi -n-m+1) 
+ 2m(5 ] ) m - l 5{yt -n-m + 2). 



") = lPE ro=2 mW kJ {Q _ n y (61) 
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SUB(l) SUB(2) SUB(3) SUB(4) 



1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 




0246824682468246 



<N> <N> <N> <N> 

Figure 1. A plot of the energy (upper panels) and the fluctuation in particle number 
(lower panels) for SUB(l) to SUB(4) ECCM based on the particle picture, for Q = 4. 
Black circles denote an exact solution. 



branch. Also, all of these solutions are infinitely degenerate; something which is not 
obvious from the diagrams! This should not surprise us: in the SUB (4) calculation each 
of these states is orthogonal to the reference state, and is (essentially) generated by a 
single A^ n excitation on the ground state. 
In general we can prove that 

(H) = \(N 2 }-^(N). (62) 

Putting in the explicit form of the ground state energy, we can show that a state where 
the expectation value of H is the ground-state energy is also a state where (N 2 ) = (N) 2 , 
i.e., AN 2 = 0. 

It should worry us that the usual implicit assumption in CCM-that the physical 
solution is smoothly connected in parameter space-is not correct in this case. If we are 
interested in a relatively low order SUB(n) calculation for a problem with many particles 
this may not be such an important problem, but this requires detailed investigation. 

Notwithstanding those concerns, if we apply the SUB (2) ECCM to a variety of 
values of Q, see Fig. [2j we see that this is quite a successful and stable improvement to 
the BCS below the half-filled shell. Above that we could work with holes relative to the 
full shell and get much better results. 

We can now concentrate on the "physical branch" of solutions, identified by the 
smooth connection to the ground state, and use solution following techniques to trace 
the higher order solutions for this value. In other word, we only solve a small subset. 
We choose SI = 10, and find the result of Fig. [31 We see that we can indeed improve the 
results considerably, but also note the indication that there is a limit to the improvement. 
This is obvious from the fact that we know this is not the exact solution as we increase 
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n=4 n=5 n=6 a=s n=io a=u n=2o 




4 6 8 2 4 6 8 10 5 10 5 10 15 5 10 15 20 5 10 15 20 25 10 20 

<N> <N> <N> <N> <N> <N> <N> 



Figure 2. A plot of the energy (upper panels) and the fluctuation in particle 
number for SUB (2) ECCM based on the particle picture, for a variety of values of 
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Figure 3. A plot AE = E — £ cxact vs (N) for il = 10 using ECCM, with "ordinary" 
operators, rather than quasiparticle ones. We compare different levels of truncation, 
with either hole operators acting on the completely filled state (full) or particle 
operators acting on the empty state. 



the truncation level. 



6.3.2. Brueckner orbitals We would like to analyse the same systems for the 
quasiparticle ECCM (and NCCM). Unfortunately getting a full set of solutions to the 
quasiparticle ECCM, even for Q = 4, as in Fig. [1] has proven to be impossible-the 
equations are just too complicated to get a complete set of solutions beyond SUB(2). 
If we look at the quasiparticle NCCM and ECCM for the SUB(2) case (only s (2) and 
s^are non-zero) we get a rather interesting looking set of solutions, see Fig. [6j 

In that figure we see the usual multitude of solutions; the physically correct one 
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Figure 4. A plot of the energy (upper panels) and the fluctuation in particle number 
(lower panel) for SUB(2) quasiparticle NCCM based on the quasiparticle picture, for 
a variety of values of the degeneracy SI. 
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Figure 5. A plot of the energy (upper panels) and the fluctuation in particle number 
for SUB(2) quasiparticle ECCM based on the quasiparticle picture, for a variety of 
values of the degeneracy SI. 



(that turns to E — for (N) = 0) collapses at mid-shell. This agrees with the fact 
that AN 2 goes negative around midshell, showing this is or has become an unphysical 
solution. Doing the same calculation with NCCM for the operators gives a better 
result, to our great surprise. 

If we increase the order of the calculation, we see a suggestion of convergence for 
the odd (but not even orders). 

We can perform a similar NCCM calculation. Here we can show that the SUB (2) 
calculation in the restricted space is exact (the full solution collapses to the "two-pair" 
form assumed in all the above calculations). The solution collapses in higher order, 
with the strangest result at order 4. The sharp turns are real-our solution following 
technique resolves these exactly. 




Figure 6. A plot E as a function of (N) for SI = 10, using the quasispin ECCM. We 
compare different levels of truncation, with either the filled (full) or empty state as a 
reference state. 
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Figure 7. A plot AE = E - E cxact vs (N) for Jl = 10 , using the quasispin ECCM. 
We compare different levels of truncation, with either the filled (full) or empty state 
as a reference state. 



6.3.3. Maximum overlap orbitals We can finally try to provide one more set of 
solutions, where we fix the SUB(l) states by a fixed mean-field calculation, and only 
allow SUB(2) excitations beyond that (if we allow further SUB(l) excitations, we are 
back to the case of Brueckner orbitals). Unfortunately, there are no solutions to the 
resulting equations, i.e., we can't both have the particle number correct at mean-field 
and at quasiparticle SUB(2) level. 

6.3-4- Comparison of methods All of this raises a couple of questions: 




Figure 8. A plot E vs n for ft = 10, using the quasispin NCCM. We compare different 
levels of truncation, with either the filled (full) or empty state (empty) as a reference 
state. 



(i) Should we use (trust) the ECCM or the NCCM based calculations-we have no clear 
answer at the moment 

(ii) Should we use the quasispin approach, or should we stick to the ordinary ECCM 
without quasiparticles? This can be answered, since the quasispin algebra make it 
easy to a calculation with 

and similar for the full shell. The results, shown in Fig. [3J, suggest that there is no 
gain in using this approach, but in general the approach using the quasiparticles is 
much more compact-especially if we stick to the NCCM truncation. 

Before drawing too strong a conclusion, we need to look at the harmonic fluctuations, 
which in this case is the solution to the generalised RPA. 

6.3.5. Harmonic excitations In order to look at the stability of the solutions and 
excitations, we must calculate the harmonic fluctuations about the equilibrium, also 
called the RPA for the Brueckner orbitals (the particle case is trivial, and was already 
discussed above). The easiest way to derive this is to use the time- dependent variational 
principle, derivable from the action 

A[s,s} = j dt(4>\e S e~ s {H - id t )e s \<j)). (63) 

The only tricky point is that is time-dependent, as are the operators appearing in S. 
The expectation value of H has already been evaluated, and we can easily find that (we 
use "NC" as a short-hand for the nested commutator expansion) 

(4>\e § e- s (td t )e s \<P) = (4>\e S e~ s {id t e s )\<t>) + (0|e^|0) 
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71-1 



+ (4>\e^isA^\(f)). 



(64) 



The first term can be trivially calculated with the techniques we have developed; the 
other two require (a bit of) further algebra. Let us start with the final term. We need 
to express A in quasispin operators, 

A f = 44 = ( Sc * + 4) + 4) =-&S + 6*- s(n - ft). (65) 



We thus find that 

is(4>\e S A^\(l>) = ft55, 

since S only contains terms with two or more <5's. 

We are left with the middle operator, which requires most work. First of all 

d t 8^ = d t |^(1 — ss) a\ — saj. (1 — ss) a| + sa k 

-a- k d\ + d\a k 



(66) 



= — (5s + ss) 
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+ s 



(ss + ss) (§c k + 4) 4 + 4 (y~ sc k + 4) 



+ 5 



(1 — ss) Cfc — sd\ 4 + 4 (l — ss)ck + sd 



= [-2(ss + ss) + 2ss] 44 
+ [— (ss + ss)s — 5(1 — ss)] 
+ [(ss + 5s)5 + 5(1- ss)] d\c k 

= - 255^+ (5 + 5 2 5) (n-ft). 
This can be used to find that 

m— 1 
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= - 255m + (5 + 5 2 5) (5^) m 1 m [n - ft + m - 1] . (68) 

This can be used to evaluate 
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Figure 9. The RPA frequencies for the normal ECCM for — 10 at various levels of 
truncation 
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Figure 10. The RPA frequencies for the quasiparticle ECCM for D, = 10 at various 
levels of truncation 
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6.3.6. Result from RPA calculations We have preformed RPA calculation for all of 
the approximation schemes discussed above; we shall concentrate on two of most 
general schemes, ECCM in normal operators and quasiparticle ECCM, Figs. M and 
IT0| respectively. These two figures should be looked at in conjunction with [3] and [3 
In each of these we concentrate on those states with the same symmetry as the ground 
state: our operators are only zero angular momentum pairs. 

From both these figures, we note that we do see the expected zero modes for 
excitations breaking particle number, but remaining on the solution branch being 
studied. These are of course doubly-degenerate. As we increase the level of truncation, 
more-and-more other modes come in. These also break particle number, as can be 
seen e.g. at (N) = 0, where the only excited states available are those with different 
particle number, e.g., the iV = 2 states for SUB (2), iV = 2 + 4 for SUB (3), etc. The 
excitations energies are of course also identical for the empty shell for the two schemes, 
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since they coincide at that point. Things change as we move away towards the middle of 
the shell. Whereas the standard ECCM gives real frequencies, the quasiparticle ECCM 
solutions meet in pairs (and turn complex), showing a break-down of the approximation 
scheme. Also, the break-down point seems to move to smaller N as we increase the 
level of truncation. It also sheds light on the rather different behaviour for even and 
odd truncation levels seen before: for an even truncation scheme, we always have one 
odd pair of modes, taking away the pair of zero modes, whereas for an odd truncation 
level all non-zero modes can and will collide in pairs. 

7. Discussion 

We have shown that one can formulate a number of CCM based methods for application 
to problems with pairing. These have been analysed for a rather simple model, and all 
seem to suffer from some problems. The particle ECCM is the most stable method for 
low order calculations, but all methods seem to suffer some drawback. The real surprise 
is that the low order particle ECCM fails as badly as it does-we would have expected 
it to work best. No detailed understanding of these results exists, but they are clearly 
linked to the restoration of particle number required in these calculations. 

Clearly the exact integrability of the underlying model may also play a role in these 
results, but one feature will persist: the restoration of the broken particle number means 
that the exact solutions have zero SUB(l) pairing, suggestive of the fact that different 
solutions lie on different branches of CCM solutions. This will also lead to instabilities 
at some order-we do expect that for real system with pairing in large model spaces the 
problems are much less pronounced, but the comparison of SUB2 NCCM and ECCM 
seems to be a promising aspect of any approach, as is the study of the particle number 
fluctuations, which seems to be the best tool to highlight unphysical solutions. 

If we were to try and apply an improved many-body method to fermions, e.g., 
confined in a harmonic trap, any one of the methods discussed here is probably a good 
choice, as long as we don't push the order of approximation too high; we must, however, 
carefully monitor the behaviour of fluctuation in particle number and the local-harmonic 
modes with ground state symmetry. 
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